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Abstract 


Steady-state  flow  and  structural  profiles  of  immiscible  components  A  and  B  (molecular 
weights  M A  and  M B,  MA<MB)  in  a  non-conservative  open  system  are  studied  by  an 
interacting  lattice  gas  Monte  Carlo  simulation.  Concentration  gradient  and  hydrostatic  bias  H 
drive  the  constituents  (A,B)  which  are  continuously  released  from  the  bottom  with  equal 
probability  against  gravity.  At  low  bias,  the  segregation  of  A  and  B  leading  to  a  partial 
layering  is  enhanced  toward  the  bottom.  The  longitudinal  density  profile  with  a  high  density  in 
the  bottom  region  and  low  toward  the  top  shows  linear,  exponential,  and  power-law  decays  in 
different  regions  of  depth  or  altitude  which  varies  systematically  with  the  pressure  bias.  The 
transverse  density  profiles  show  segregation  with  different  domain  sizes  and  layering 
depending  on  the  bias.  Response  of  their  steady-state  flux  density  j  to  the  hydrostatic  bias 
H  is  found  to  be  linear  at  higher  bias.  The  difference  in  response  of  the  flux  density  of  the  two 
components  becomes  more  pronounced  at  low  bias  and  higher  miscibility  gap. 

Published  by  Elsevier  B.V. 

Keywords:  Lattice  gas;  Mixture;  Segregation;  Non-equilibrium  steady  state;  Linear  response 


A  continuous  release  of  constituents  (particles  (A,B))  from  a  source  into  the 
bottom  of  an  effective  open  medium  leads  to  its  steady-state  flow.  How  much  the  net 
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amount  of  A  and  B  flows  and  how  their  density  profiles  appear,  depends  not  only  on 
the  rate  of  release  but  also  on  the  type  of  constituents  (i.e.,  immiscibility),  their 
interaction  with  the  effective  medium  (i.e.,  miscibility),  sedimentation  and 
precipitation  (due  to  molecular  weight),  pressure  bias  (driving  field),  etc.  In  order 
to  gain  a  systematic  understanding  of  evolving  morphologies  and  their  flow,  one  has 
to  narrow  down  the  parameter  space.  For  example,  we  can  fix  the  rate  of  release  of 
constituents  from  the  source  and  their  molecular  weights  (i.e.,  MA/MB  =  ^)  and  vary 
the  pressure  bias  ( H ).  One  may  study  the  effect  of  miscibility  by  changing  the 
interaction  strength  among  the  constituents  (AA,AB,BB)  and  with  the  effective 
medium  (AS,BS),  i.e.,  with  the  empty  lattice  sites  ( S ).  Identifying  the  qualitative 
trend  from  the  unsealed  quantitative  data  generated  by  computer  simulations  is  one 
of  our  objectives  in  this  article. 

The  global  (and  local)  structure  evolving  from  individual  interacting  components 
and  their  flow  share  some  general  characteristics,  at  least  in  part  with  many  systems 
such  as  flow  of  immiscible  fluid  mixture  (oil  and  water)  into  a  miscible  solvent,  flow 
of  methane  gas  and  hydrocarbon,  and  dissociation  of  hydrates  in  aqueous 
environments  [1].  Our  main  interest,  however,  is  focused  on  understanding  the  flow 
of  a  mixture  of  methane  gas  along  with  other  components  (i.e.,  immiscible 
hydrocarbon)  below  the  ocean  floor  to  sub-bottom  area  in  aqueous  and  dynamic 
sedimentary  environment  [2-5].  In  the  long  term,  we  are  interested  in  a 
comprehensive  understanding  of  the  effects  of  driving  fields  such  as  concentration 
gradients,  hydrostatic  pressure,  thermal  gradient,  etc.  [6,7]  which  are  crucial  in 
understanding  the  hydrate  formation,  mound  evolution,  dissociation  leading  to  such 
instabilities  as  mud  volcanoes  [8-10]  below  the  ocean  floor.  We  recently  reported  [11] 
the  onset  of  segregation  with  partial  layered  patterns  in  a  steady-state  flow.  Changes 
in  concentration/density  profiles  with  the  hydrostatic  pressure  bias  were  also 
examined  very  recently  [12]  for  immiscible  mixtures  (A,  B)  with  a  range  of  molecular 
weight  ratios  (MB/MA  =  1-10).  Attempts  were  made  to  quantify  the  variations  of 
the  overall  steady-state  density  of  these  constituents  and  their  scaling  with  the  bias 
(//)  and  the  molecular  weight.  In  this  article,  we  study  the  response  of  flow  to  bias 
with  steady-state  structures  involving  segregation  where  both  constituents  are 
released  from  the  source  regardless  of  their  overall  concentration  in  the  system. 

The  mixture  consists  of  interacting  particles  which  are  driven  by  hydrostatic  bias 
and  concentration  gradient  against  gravitational  forces  while  they  aggregate  and 
segregate.  Such  non-equilibrium  steady-state  systems  provide  a  wide  range  of 
structural  patterns.  Effects  of  external  driving  field/force  on  ordering  and  non¬ 
equilibrium  phase  transition  in  conservative  particles  systems  (with  constant  number 
of  constituents)  are  studied  extensively  with  a  long  list  of  literature  [13,14]. 
Investigations  of  structures  and  patterns  in  granular  driven  systems  [15,16]  have 
attracted  a  considerable  interests  in  recent  years  [17-23].  Understanding  the  global 
patterns  from  the  mobility  of  individual  constituents  are  difficult  to  monitor 
experimentally  [24].  However,  it  is  relatively  easier  to  investigate  flow  and  patterns 
[25,26]  in  such  a  complex  particulate  mixture  by  computer  simulation  methods  such 
as  lattice  gas  [27,28],  molecular  dynamics  [29,30],  and  Monte  Carlo  [31]  methods. 
Most  of  the  particulate  driven  system  [13-23]  deal  with  conservative  systems  with  a 
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constant  number  of  constituents.  Unlike  Boltzmann  lattice  gas  [32],  it  is  easier  to 
incorporate  interactions  between  constituent  particles  and  effective  medium  (empty 
lattice  sites,  see  below)  in  interacting  lattice  gas.  We  use  an  interacting  lattice  gas 
Monte  Carlo  simulation  here  to  study  the  flow  response  in  a  steady-state  mixture 
with  non-conservative  constituents  as  follows. 

As  before  [1 1,26],  we  consider  a  cubic  lattice  of  size  l)  with  L  =  30-400  with  a 
source  of  particles  A  and  B  connected  to  the  bottom  of  the  lattice  (z  =  1).  About  half 
of  the  lattice  sites  are  initially  occupied  by  A  or  B  randomly  with  equal  probability. 
The  empty  sites  represent  the  effective  solvent  ( S )  medium.  The  steady-state  global 
properties  do  not  depend  on  the  initial  configuration  (see  below).  The  interaction 
energy  U  between  these  components  is  described  by 

-w.*)  ’  o) 

i  k 

where  index  i  runs  over  all  particles  and  k  over  all  nearest  neighbor  sites  of  i  with, 
J(A,A)  =  J(B,B)  =  —J(A,  B)  =  —J(B,A)  =  -e{  ,  (2) 

J{A,S)  =  J{B,S)  =  -e2.  (3) 

The  positive  miscibility  gaps  (e\  =  1,2)  considered  here,  enhance  the  phase 
separation  between  A  and  B  in  contrast  to  their  mixing  with  the  solvent  £2  =  1. 

Downward  fall,  i.e.,  sedimentation  is  controlled  by  the  molecular  weights  of 
constituents  (M A  =  0.1,  MB  =  0.3)  via  change  in  the  gravitational  potential  energy 
as  they  move  down  (— z)  or  up  (-fz).  In  addition  to  high  concentration  gradient  at  the 
bottom  due  to  source,  a  hydrostatic  pressure  bias  (//)  is  considered  to  drive  the 
constituents  upward  against  gravity.  The  Metropolis  algorithm  is  used  to  move 
randomly  selected  particles  to  their  nearest  neighbor  sites  selected  randomly  with 
following  probabilities  along  the  transverse  (±x,  ±y)  and  vertical  (±z)  directions: 

PX  =  P_X  =  P  =  />=!;  p,  =  1  +  — ,  p_:  =  }—Ji;o^H^\.  (4) 

6  6  6 

Periodic  boundary  conditions  are  used  along  the  transverse  direction;  both  top 
(z  =  L)  and  bottom  (z  =  1)  of  the  sample  are  open  for  particles  to  escape.  A  particle 
(A  or  B  with  equal  probability)  is  released  into  a  bottom  lattice  site  as  soon  as  it  is 
vacated  by  the  particle  when  it  moves  up  or  diffuses  laterally.  As  usual,  the  attempt 
to  move  each  particle  once  defines  the  unit  Monte  Carlo  step  (MCS)  time. 

As  the  simulation  proceeds,  particles  (A,B)  move  in  and  out  of  the  sample.  The 
number  of  particles  and  their  concentration  profiles  (planar  densities  along  x ,  y,  z) 
change  but  eventually  reach  a  steady-state — a  non-equilibrium  steady-state  system. 
Simulations  are  performed  with  many  independent  samples  each  for  a  sufficiently 
long  time  steps  to  reach  steady  state.  We  have  observed  [11]  evolution  of  interesting 
structural  patterns,  i.e.,  onset  of  phase  separation,  oscillation  in  lateral  density 
profiles,  layering,  etc.  These  structural  patterns  share  some  characteristics  of 
granular  systems  studied  in  recent  years  [33-35].  In  the  following,  we  focus  mainly  on 
the  global  transport  and  flow  response  in  segregating  steady-state  patterns. 
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Simulations  are  performed  on  various  lattice  sizes  303-4003  with  most  production 
runs  on  303,  503,  and  1003  lattices  using  as  many  as  128  independent  runs  for  the 
range  of  pressure  bias  H  =  0.0- 1.0.  Finite  size  effects  on  the  qualitative  behavior  of 
the  physical  quantities  are  negligible  particularly  with  503  and  1003  samples  except 
for  few  runs  at  low  bias  values.  The  miscibility  gaps  are  restricted  to  e\  =  1.0,  2.0, 
and  £2  =  1.0  with  the  molecular  weights  MA  =0.1,  Mb  =  0.3,  and  temperature 
T  =  1.0  in  arbitrary  units  of  interaction  energy  and  Boltzmann  constant. 

Fig.  1  shows  the  steady-state  density  profiles  in  the  whole  range  of  pressure  bias. 
The  longitudinal  profiles  with  altitude  clearly  show  a  decay  of  the  density  with  the 
height  (z):  sharp  decline  at  low  height  is  followed  by  a  slow  but  linear  decay  at  the 
intermediate  height  before  a  sharp  decline  again  toward  the  top.  The  density  profiles 
depend  on  the  magnitude  of  the  bias.  The  density  of  heavier  component  ( B )  is 
generally  higher  particularly  at  lower  bias  and  low  altitudes.  At  extreme  values  of 
bias  (H  ->  1),  the  density  profiles  of  both  components  approach  the  same  shape. 
Attempts  are  made  to  explore  the  exponential  (log-normal  scale)  and  log-log  (Fig.  2) 
fits  which  suggest  that  the  steady-state  variation  of  the  density  of  each  component 
depend  on  the  range  of  altitude,  magnitude  of  the  bias,  and  the  molecular  weight  of 
the  components.  For  example,  at  the  high  value  of  the  bias  (//  =  1.0)  where  the 
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Fig.  1.  Density  ( dA,dB )  profiles  of  A  and  B  components  at  a  range  of  bias  H  =  0.15-1.00.  The  inset 
shows  a  transverse  (F)  density  profile  at  H  =  0.15.  Sample  size  1003  is  used  with  32-128  runs. 
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Fig.  2.  Density  (< dA,dB )  profiles  of  A  and  B  components  at  a  range  of  bias  H  =  0.15-1.00  on  log-log  scale 
to  search  for  a  power-law  dependence.  Sample  size  1003  is  used  with  32-128  runs  (Fig.  1  on  a  different 
scale). 


density  profiles  of  A  and  B  are  similar  (see  Fig.  2),  a  sharp  power-law  decay  of 
density  with  the  altitude  (at  Z2:  1—8)  is  followed  by  a  slower  power-law  decay 
(Z21  15-80)  before  dropping  precipitously  toward  the  top  (Z  2:  90- 100).  The 
power-law  decay  (Fig.  2)  at  H  =  1.0,  dA/B  oc  z-04  at  lower  altitude  is  followed  by 
dA/B  oc  z~ 02  at  higher  altitude.  Note  that  the  power-law  exponent  depends  on  the 
pressure  bias  and  range  of  the  altitude.  Particularly  at  low  bias  values,  a  range  of 
power-law  decays  including  non-monotonic  variation  of  densities  with  z  is  seen 
where  concentration  gradient  seems  to  dominate  the  driving  forces. 

Toward  the  bottom  end  where  both  components  have  higher  densities,  it  is 
easier  to  identify  the  structural  changes  due  to  their  interaction  (i.e.,  miscibility). 
The  density  profile  is  important  as  its  gradient  provides  one  of  the  driving  bias 
for  the  mass  transfer  from  bottom  to  top.  The  transverse  density  profile  provides 
insight  into  the  lateral  distribution  of  their  constituents  which  affect  their  transverse 
mobility  (see  below).  In  fact,  the  phase  separation  between  A  and  B  is  most  visible 
at  low  values  of  bias  (H  =  0.0-0.30);  a  transverse  profile  is  shown  in  the  inset  of 
Fig.  1  at  H  =  0.15.  The  complementary  and  somewhat  oscillatory  variation  in  the 
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density  of  B  (higher)  and  A  (lower)  exhibits  segregation  with  onset  of  their  layering 
toward  the  bottom  planes  (also  seen  in  visualization).  The  segregated  layers  smear 
out  toward  the  upper  half  where  their  density  become  much  smaller. 

Although  each  particle  (A,B)  attempts  to  move,  particles  at  the  top  (low-density 
phase)  are  much  more  mobile  than  those  toward  the  bottom  (high-density  phase). 
The  motion  of  each  particle  may  not  be  as  insightful  as  the  overall  motion  of  their 
center  of  mass  for  the  net  flow.  Fig.  3  shows  variation  of  the  longitudinal  component 
of  the  mean  square  displacement  on  a  log-log  scale  for  the  whole  range  of  bias; 
variation  of  the  representative  transverse  (y)  component  is  presented  in  the  inset.  The 
power-law  behavior  is  observed  even  at  the  lowest  value  the  bias  (//  =  0.15)  in  the 
asymptotic  time  limit.  Despite  similarity  in  the  power-law  dependence,  there  are 
considerable  differences  in  their  magnitude  (R2);  i.e.,  it  is  higher  for  the  lighter  (A) 
component  particularly  at  low  bias.  However,  the  variation  with  bias  is  more 
pronounced  for  the  heavier  component  (B).  Note  the  contrast  in  the  order  of 
magnitude  of  longitudinal  (R2)  and  transverse  (R2)  components:  the  extremely  low 
values  of  R2y  show  very  low  overall  mobility  of  the  center  of  mass  along  the 
transverse  direction  while  the  comparatively  high  magnitude  of  R2  show  a  large 
transport  along  vertical  direction. 


Fig.  3.  Mean  square  displacement  (longitudinal  component  /??)  versus  time  steps  for  components  A  and 
B\  variation  of  a  transverse  (R2V)  component  is  presented  in  the  inset  at  extreme  values  of  H  =  0.15, 1.00. 
Statistics  is  the  same  as  Fig.  1. 
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The  constant  in-flux  of  constituents  from  the  source  at  the  bottom  in  presence  of 
an  upward  pressure  bias  provides  their  net  flow  (i.e.,  mass  transfer)  despite 
gravitational  forces.  In  steady  state,  the  output  from  the  top  must  be  equal  to  net 
input  flow  at  the  bottom  as  the  conservation  of  mass  leads  to  continuity  equation. 
The  rate  (j)  of  mass  transfer  per  unit  cross-sectional  area  ( Ac )  defines  the  flux  density, 


J_dQ 
Ac  d t 


(5) 


where  Ac  =  Lx  x  Ly  and  Q  is  the  net  mass  transfer  at  time  step  t.  In  the  steady  state 
the  flux  density  from  the  top  and  bottom  should  be  equal.  Variation  of  the  flux 
density  with  the  time  steps  for  both  components  is  presented  in  Fig.  4  which  is 
consistent  with 


Jtop  J  bottom  (b) 

in  the  long  time  limit.  The  approach  to  steady  state  depends  on  the  bias,  the  larger  is 
the  bias,  the  faster  is  the  approach  to  steady  state. 

One  would  expect  the  flux  density  to  increase  on  raising  the  driving  bias.  In 
addition  to  the  pressure  bias  (//),  the  concentration  gradient  also  drives  the 


jA 


Fig.  4.  Flux  density  versus  time  steps  at  various  bias.  Upper  data  points  are  the  flux  density  from  the 
bottom  and  lower  points  are  for  flux  density  from  the  top.  Statistics  is  the  same  as  Fig.  1 . 
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constituents  upward  against  gravitational  force.  Since,  the  downward  gravitational 
bias  (which  depends  on  the  molecular  weight  of  constituents)  is  constant,  and  it  is 
difficult  to  quantify  the  overall  bias  due  to  concentration  gradient,  the  response  of 
the  flux  density  j  to  the  driving  field  is  examined  by  varying  the  pressure  bias  (//). 
Fig.  5  shows  the  variation  of  the  flux  density  j  with  bias  //,  for  different  sample  sizes 
and  interaction  parameter  s\  =  1,2  for  the  miscibility  gap.  In  the  high  range  of  bias, 
the  response  of  j  to  H  is  linear  and  is  less  sensitive  to  miscibility  gap.  In  the  low-bias 
regime,  there  is  deviation  on  increasing  the  miscibility  gap  from  e\  =  1  to  2.  The 
response  at  low  values  of  H  differs  from  that  in  the  high  range  of  bias.  Effects  of 
gravitational  forces  and  the  concentration  gradient  are  clearly  more  visible  at  low 
values  of  pressure  bias.  Further,  the  effect  of  molecular  weight  ( Ma<Mb )  on  the 


Fig.  5.  Variation  of  the  steady-state  flux  density  with  the  pressure  bias  H  with  different  sample  sizes 
(L  =  30,50, 100)  and  miscibility  gap  ei  =  i—  1,2.  The  number  of  independent  runs  varies  from  10  (L  = 
30,50)  to  128  (L=  100)  sizes. 
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response  is  also  significant  particularly  at  low-pressure  bias:  A  component  responds 
faster  than  B  as  jA  ^ jB .  Within  the  statistical  fluctuations,  there  is  hardly  any  finite 
size  effect  on  the  qualitative  nature  of  the  response. 

In  summary,  an  interacting  lattice  gas  computer  simulation  model  is  used  to  study 
the  flow  response  of  immiscible  components  ( A,B )  which  are  constantly  released 
with  equal  probability  from  the  source  at  the  bottom  and  driven  by  an  upward  bias 
against  gravitational  force.  The  number  of  particles  is  not  conserved  as  they  escape 
from  the  top  and  leak  from  the  bottom.  These  interacting  constituents  in  an  effective 
medium  are  continuously  driven  out  of  equilibrium  by  the  bias  but  their 
concentration  profile  and  flux  density  reach  steady  state.  The  longitudinal  density 
profile  with  a  high  density  in  the  bottom  region  and  low  toward  the  top  shows  linear, 
exponential,  and  power-law  decays  in  different  regions  of  depth  or  altitude  which 
varies  systematically  with  the  pressure  bias.  The  transverse  density  profiles  show 
segregation  with  different  domain  sizes  depending  on  the  bias.  The  flux  density 
( JaJb )  of  both  constituents  respond  linearly  to  pressure  bias  in  the  high-bias  regime 
in  which  jA^jB  as  H  ->  1.  The  difference  in  flow  response  increases  on  reducing  the 
bias.  In  the  low-bias  region,  the  response  of  flux  density  of  the  component  ( B )  with 
the  higher  molecular  weight  is  lower  than  that  of  the  lighter  constituents  (A).  The 
flux  density  response  is  also  found  to  be  lower  with  higher  miscibility  gap  in  the  low- 
bias  region  where  not  only  the  response  is  different  from  that  in  the  high-bias  region, 
but  also  signals  a  deviation  from  a  linear  response. 
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